subroutine srccsd(ie, fij, fab, fia, & 
          iijkl, iijka, iabij, &
          iiabj, iabci, iabcd, &
          t1, t2, r1, r2, &
          fae, fmi, fme, &
          wmnij, wmbej, icore, & 
          maxcore, escf, eccsd) 
use mod_iop
use mod_size
use mod_orbit
use mod_ioff
implicit none
include 'omp_lib.h'
!
integer prec2, isd2, igpu2, totgpu2
common/precsoc/prec2, isd2, igpu2, totgpu2
!
integer, intent(in) :: maxcore 
real*8, intent(inout) :: escf, eccsd
real*8, intent(inout) :: ie(*), fij(*), fab(*), fia(*) 
real*8, intent(inout) :: t1(*), t2(*), r1(*), r2(*) 
real*8, intent(inout) :: iijka(*), iabij(*), iijkl(*)
real*8, intent(inout) :: fae(*), fmi(*), fme(*)
real*8, intent(inout) :: wmnij(*), wmbej(*)
real*8, intent(inout) :: iiabj(*), iabci(*), iabcd(*), icore(*)
!
real*4, allocatable,device :: icoreg(:)
!
integer nsize, i, j, a, b, icrsizg
integer maxvar, maxdim, iter
integer i0, i1, i2, i3, i4, i5, i6, i7, i8, i9
integer i10, i11, i12, i13
real*8 dtmp, errdiis, eccsd2, t1diag
real*8 eso, tt(20), rt(20), tt1, tt2
real*8 Eccsdio(1000), iter_time(1000)
real*4, allocatable :: idoso(:), itmp(:)
!
real*8, external :: ddot, dnrm2, nrm2so
integer, external :: idsymsoc, irpdso
!
zorb(2:8) = 0
if(totgpu==1) then 
   icrsizg = 1024*1024*1024/4*11 ! 11GB
   icrsizg = igpusize ! 11GB
   allocate(icoreg(icrsizg))
endif
!ldebug = .false.
!ldebug = .true.
!
if(ldebug) then 
  write(6,'(A,8I4)')'vrta', vrta
  write(6,'(A,8I4)')'vrtb', vrtb
  write(6,'(A,8I4)')'popa', popa
  write(6,'(A,8I4)')'popb', popb
  write(6,*)'nioi3', nioi3
  write(6,*)'nioi3', idsymsoc(1, vrta, vrta, 0, popa, popa, 0)  
  write(6,*)'nioi3',   idsymsoc(1, vrta, vrtb, 1, popa, popb, 1)   
  write(6,*)'nioi3',  idsymsoc(1, vrtb, vrtb, 0, popb, popb, 0)   
endif 
call myicopyso(iabij, r2, nioi3)
call myicopyso(fia, r1, nifvo)

if(ldebug) then 
call printdnrm2('iabij1', iabij, 3)
call printdnrm2('iabij1', r2, 3)
endif
!
maxvar = nifvo + nioi3
if(ioiabcd==0) then 
  allocate(idoso(maxvar*12*isd), itmp(maxvar*isd)) 
endif
iter = 0 
if(ldebug) then 
  call printdnrm2('iabij', iabij, 3)
endif
!
!--iteration begins
!
do
iter = iter+1
write(6,*)
write(6,*)'=============in iteration:', iter
!
call cpu_time(tt(1)); call walltimeso(rt(1))
call denomt1(t1, r1, ie, 1)
call denomt2(t2, r2, ie, 1, 1, 1)
if(ldebug) then 
  call printdnrm2('r2 in the beginning', r2, 3)
  call printdnrm2('t2 in the beginning', t2, 3)
endif
!
maxdim = 5
call cpu_time(tt1)
if(ioiabcd==0) then 
call dosodiis3(iter, t1, maxvar, maxdim, icore, &
              maxcore+maxvar, errdiis, idoso, itmp)
else
  call dosodiis2(iter, t1, maxvar, maxdim, icore, &
                maxcore+maxvar, errdiis)
endif
call cpu_time(tt2)

if(ldebug) then 
 write(6,*)'time for dosodiis2', tt2-tt1
endif
call effectt2(t2, t1,  icore, 1.d0, no, nv)
if(ldebug) then 
call printdnrm2('t2after doso ', t2, 3)
endif
!
call calcene(t1, t2, fia, iabij, eccsd)
call calcene2(t1, t2, fia, iabij, eccsd2)
eccsdio(iter)=eccsd2
write(6,*)'CCSD energy', eccsd2
!write(6,*)'CCSD energy', eccsd, eccsd2
!write(6,*)'total energy', escf+eccsd, escf+eccsd2
!
call effectt2(t2, t1, icore, -1.d0, no, nv)
if(iter>maxiter) exit
!
call cpu_time(tt(2)); call walltimeso(rt(2))
if(ldebug) then 
 write(6,'(A, F10.3, F10.3)')'====time for energy     ', tt(2)-tt(1), rt(2)-rt(1)
endif
!===========prapere for the begin 
call myicopyso(fab, fae, nifvv)
call myicopyso(fij, fmi, nifoo)
call myicopyso(fia, fme, nifvo)
call myicopyso(iijkl, wmnij, nioi1)
call myicopyso(iiabj, wmbej, nioi4)
call myicopyso(fia, r1, nifvo)
call myicopyso(iabij, r2, nioi3)
call cpu_time(tt(3)); call walltimeso(rt(3))
if(ldebug) then 
   write(6,'(A, F10.3, F10.3)')'====time for prepare    ', tt(3)-tt(2), rt(3)-rt(2)
endif
   !
!==============n5 term=====================
!
if(ldebug) then 
 call printdnrm2('iiabj', iiabj, 4)
 call printdnrm2('t1-start' , t1, 22)
 call printdnrm2('r1-start' , r1, 22)
endif
call calciiabj(r1, r2,  t1, &
          iiabj, icore)
 call cpu_time(tt(4)); call walltimeso(rt(4))
if(ldebug) then 
  call printdnrm2('r1-iiabj' , r1, 22)
  write(6,'(A, F10.3, F10.3)')'====time for n5 term    ', tt(4)-tt(3), rt(4)-rt(3)
endif
!=============iabci term==========
!call printdnrm2('iabci' , iabci, 5)
!ioiabci=1
if(ioiabci==1) then 
 call calciabci_n5(r1,r2,t1,t2,fae,wmbej,icore)
elseif(ioiabci==0) then 
 call r1ai_t2efmi_iefam(r1, t2, iabci, icore, -0.5d0)
 call r2abij_t1ei_iabej(r2, t1, iabci, icore, 1.d0)
 call wmbej_t1fj_iefbm(wmbej, iabci, &
                   t1, fae, icore, -1.d0, 1.d0) 
endif
igpu=0; igpu2=0
call cpu_time(tt(5)); call walltimeso(rt(5))
if(ldebug) then 
  write(6,'(A, F10.3, F10.3)')'====time for wmbej large', tt(5)-tt(4), rt(5)-rt(4)
  call printdnrm2('r1-t2abci' , r1, 22)
  call printdnrm2('r2-n5term' , r2, 3)
  call printdnrm2('fae', fae, 11)
  call printdnrm2('wmbej', wmbej, 4)
endif
!==============================
!-----terms which are independent 
!=============================
!
!==========abci======================
call cpu_time(tt(6)); call walltimeso(rt(6))
if(totgpu==1) then
  i0=1
  i1 = i0 + nifvo*isd
  i2 = i1 + nioi3*isd
  i3 = i2 + nifvo*isd
  i4 = i3 + nioi3*isd
  i5 = i4 + nvvoo*isd*3
  igpuleft = igpusize - i5
  if(i5>icrsizg) goto 200
  igpu=1; igpu2=1
  call myicopy_g(r2,    icoreg(i1), nioi3)
  call myicopy_g(t1,    icoreg(i2), nifvo)
  call myicopy_g(t2,    icoreg(i3), nioi3)
  write(6,*)'calc iabci large in GPU'
  call calciabcig(icoreg(i1), icoreg(i2), & 
                 icoreg(i3), icoreg(i4), iabci, icore)
  call myicopy_c(icoreg(i1), r2,    nioi3)
else
200 if(totgpu==1) then 
     write(6,'(A, I25, I15, A)')'gpu memory not enough for iabci ', & 
                               i6, i6/1024/1024*4, 'MB'
    endif
   call calciabci(r1, r2, t1, t2, iabci, icore)
endif
call cpu_time(tt(7)); call walltimeso(rt(7))
write(6,'(A, F10.3, F10.3)')'====time for iabci      ', tt(7)-tt(6), rt(7)-rt(6)
  igpu=0; igpu2=0
if(ldebug) then 
  call printdnrm2('r1-iabci' , r1, 22)
  call printdnrm2('r2-iabci', r2, 3)
endif
!
!==========abcd======================
if(totgpu==1) then
  i0=1
  i1 = i0 + nifvo*isd
  i2 = i1 + nioi3*isd
  i3 = i2 + nifvo*isd
  i4 = i3 + nioi3*isd
  i5 = i4 + nvvoo*isd
  if(i5>icrsizg) goto 300
  igpuleft = igpusize-i5
  igpu=1; igpu2=1
  write(6,*)'calc iabcd large in GPU'
  call myicopy_g(r1,    icoreg(i0), nifvo)
  call myicopy_g(r2,    icoreg(i1), nioi3)
  call myicopy_g(t1,    icoreg(i2), nifvo)
  call myicopy_g(t2,    icoreg(i3), nioi3)
!  write(6,*)'finished copy of iabcd'
  call calciabcd(icoreg(i0), icoreg(i1), icoreg(i2), & 
                 icoreg(i3), iabcd, icoreg(i4))
  call myicopy_c(icoreg(i0), r1,    nifvo)
  call myicopy_c(icoreg(i1), r2,    nioi3)
!  call calciabcd(r1, r2, t1, t2, iabcd, icore)
elseif(totgpu==0) then 
300 if(totgpu==1) then 
     write(6,'(A, I25, I15, A)')'gpu memory not enough for iabcd ', & 
                               i5, i5/1024/1024*4, 'MB'
    endif
  call calciabcd(r1, r2, t1, t2, iabcd, icore)
endif
igpu=0; igpu2=0
call cpu_time(tt(8)); call walltimeso(rt(8))
if(ldebug) then 
  call printdnrm2('r1-iabcd', r1, 22)
  call printdnrm2('r2-iabcd', r2, 3)
endif
write(6,'(A, F10.3, F10.3)')'====time for iabcd      ', tt(8)-tt(7), rt(8)-rt(7)
call cpu_time(tt(9)); call walltimeso(rt(9))
if(ldebug) then 
 call printdnrm2('r1', r1, 22)
 call printdnrm2('r2', r2, 3)
 call printdnrm2('r1_b_wmbej' , r1, 22)
endif
!
!=======    wmbej   ========================================
!
if(totgpu==1) then 
  i0=1                  !r1
  i1 = i0 + nifvo*isd   !r2
  i2 = i1 + nioi3*isd   !t1
  i3 = i2 + nifvo*isd   !t2 
  i4 = i3 + nioi3*isd   !iijka
  i5 = i4 + nioi2*isd   !iabij
  i6 = i5 + nioi3*isd   !wmbej
  i7 = i6 + nvvoo*isd   !icore 
  i8 = i7 + nvvoo*isd*3   
!  goto 400 
  if(i6+nvvoo*isd*4>icrsizg) goto 400
    write(6,*)'calc wmbej 2 in GPU'
if(ldebug) then 
    call printdnrm2('r2', r2, 3)
    call printdnrm2('t2', t2, 3)
    call printdnrm2('iabij', iabij, 3)
endif
    igpu=1; igpu2=1
    call myicopy_g(r1,    icoreg(i0), nifvo)
    call myicopy_g(r2,    icoreg(i1), nioi3)
    call myicopy_g(t1,    icoreg(i2), nifvo)
    call myicopy_g(t2,    icoreg(i3), nioi3)
    call myicopy_g(iijka, icoreg(i4), nioi2)
    call myicopy_g(iabij, icoreg(i5), nioi3)
    call calcwmbejg(icoreg(i4), icoreg(i5), &
          icoreg(i2), icoreg(i3), icoreg(i0), icoreg(i1), &
          icoreg(i6), icoreg(i7), wmbej, icore)
    call myicopy_c(icoreg(i0), r1,    nifvo)
    call myicopy_c(icoreg(i1), r2,    nioi3)
else 
400 if(totgpu==1) then 
     write(6,'(A, I125, I15, A)')'gpu memory not enough for wmbej 2 ', & 
                               i8, i8/1024/1024*4, 'MB'
    endif
    if(ldebug) then 
       call printdnrm2('iijka', iijka, 2)
    endif
    call calcwmbej(iijka, iabij, &
                t1, t2, r1, r2, &
                wmbej, icore)
endif
  igpu=0; igpu2=0
call cpu_time(tt(10)); call walltimeso(rt(10))
if(ldebug) then 
  call printdnrm2('r1' , r1, 22)
  call printdnrm2('r2-wmbej', r2, 3)
endif
write(6,'(A, F10.3, F10.3)')'====time for wmbej      ', tt(10)-tt(9), rt(10)-rt(9)
!==============left============================

if(totgpu==1) then 
  i0=1                  !r1
  i1 = i0 + nifvo*isd   !r2
  i2 = i1 + nioi3*isd   !t1
  i3 = i2 + nifvo*isd   !t2 
  i4 = i3 + nioi3*isd   !iijka
  i5 = i4 + nioi2*isd   !iabij
  i6 = i5 + nioi3*isd   !iiabj
  i7 = i6 
  i8 = i7 + nioi1*isd   !fae 
  i9 = i8 + nifvv*isd   !fmi
  i10= i9 + nifoo*isd   !fme
  i11= i10+ nifvo*isd   !icore
  i12= i11+ nifvo*isd
  i13= i12+ nvvoo*3*isd
  if(i13>icrsizg) goto 500
  if(ldebug) then 
     write(6,*)'calc leftterm in GPU'
  endif 
  igpu=1; igpu2=1
  call myicopy_g(r1,    icoreg(i0), nifvo)
  call myicopy_g(r2,    icoreg(i1), nioi3)
  call myicopy_g(t1,    icoreg(i2), nifvo)
  call myicopy_g(t2,    icoreg(i3), nioi3)
  call myicopy_g(iijka, icoreg(i4), nioi2)
  call myicopy_g(iabij, icoreg(i5), nioi3)
  call myicopy_g(wmnij, icoreg(i7), nioi1)
  call myicopy_g(fae, icoreg(i8), nifvv)
  call myicopy_g(fmi, icoreg(i9), nifoo)
  call myicopy_g(fme, icoreg(i10), nifvo)
  call myicopy_g(fia, icoreg(i11), nifvo)
  call calcsmall(icoreg(i4), icoreg(i5), &
            icoreg(i2), icoreg(i3), icoreg(i0), icoreg(i1), &
            icoreg(i8), icoreg(i9), icoreg(i10), icoreg(i11),&
            icoreg(i7), icoreg(i12))
  call myicopy_c(icoreg(i0), r1,    nifvo)
  call myicopy_c(icoreg(i1), r2,    nioi3)
else
500 if(totgpu==1) then 
     write(6,'(A, I125, I15, A)')'gpu memory not enough for left term ', &
                               i8, i8/1024/1024*4, 'MB'
    endif
 
call calcsmall(iijka, iabij, &
          t1, t2, r1, r2, &
          fae, fmi, fme, fia,&
          wmnij, icore)
endif
igpu=0; igpu2=0
call cpu_time(tt(11)); call walltimeso(rt(11))
  if(ldebug) then 
     write(6,'(A, F10.3, F10.3)')'====time for left       ', tt(11)-tt(10), rt(11)-rt(10)
 endif
!
!----
if(ldebug) then 
call printdnrm2('r1 end', r1, 22)
call printdnrm2('r2 end', r2, 3)
endif
!
write(6,'(A, F10.3, F10.3)')'====time for iteration  ', tt(11)-tt(1), rt(11)-rt(1)
iter_time(iter) = rt(11)-rt(1)
!
enddo
write(6,*)'the final energy is', Eccsdio(iter)
open(555,file='CCSD_ENERGY')
write(555,*) ESCF, Eccsdio(iter)
close(555)
!
!do i = 1, iter
!   write(6,'(A, I4, 2F18.12)')'ECCSD with ITER', i, Eccsdio(i), iter_time(i)
!enddo

      WRITE(6,99)
99    FORMAT(T3,'     Summary of iterative solution of CC equations ')
      WRITE(6,100)
100   FORMAT(6X,59('-'))
      WRITE(6,101)
101   FORMAT(T24,'Correlation',T46,'Total',/,T8,'Iteration',T26, &
             'Energy',T46,'Energy')
      WRITE(6,100)
      DO I=1,ITER
       WRITE(6,1000) I-1, Eccsdio(I), Eccsdio(I)+ESCF, ' DIIS '
1000   FORMAT(T10,I4,T19,F18.12,T39,F19.12,T59,A6)
      ENDDO
      WRITE(6,100)

      T1DIAG = 0.D0
      T1DIAG = nrm2so(irpdso(1,nva, noa, 1), T1, 1)**2
      T1DIAG = DSQRT(T1DIAG/dble(NOA))
      WRITE(6,*) '     T1 Diagnostic =', T1DIAG
!
      WRITE(6,103) ESCF
      WRITE(6,102) ECCSD + ESCF
103   FORMAT('        E(SCF)            = ',F20.12)
102   FORMAT('        E(CCSD)           = ',F20.12)

!
return
end 
